function helperGetCovariance(hDp, prm)
    # Create spatial covariance matrices Rg for all groups from hDp
    # Inputs: hDp[] - channel matrices per user
    #         prm   - parameter dict
    # The spatial covariance statistics assumes ergodicity in space and
    # frequency, so this will accumulate over all receiver antennas and over
    # all subcarriers. There is no temporal averaging.

    G = prm["numGroups"]
    Rg = Array{Any}(undef, G)
    Nt = size(hDp[1], 2)  # Number of transmit antennas

    # Get spatial covariance for each group
    for g in 1:G
        R = zeros(Nt, Nt)  # Define the covariance matrix R

        # Loop through all users within the group and calculate R from the
        # channel estimates of all the antennas across all subcarriers for all
        # users in group
        usersInGroup = findall(x -> x == g, prm["groups"])
        for u in usersInGroup
            Hu = hDp[u]  # Get channel estimates for this user
            Nrx = size(Hu, 3)  # Number of receive antennas
            activeCarriers = findall(row -> all(!isnan, Hu[row, :, 1]), 1:size(Hu, 1))

            # Accumulate estimates over all antennas for this user
            for r in 1:Nrx
                # Accumulate estimates over all subcarriers
                for sc in activeCarriers
                    h = collect(transpose(Hu[sc, :, r]))
                    R += h * h'
                end
            end
        end
        Rg[g] = R / maximum(R)  # Normalize to get good rank measurement
    end

    return Rg
end

# [EOF]